/*
 * AnalyzeCharges.h
 *
 *  Created on: May 26, 2011
 *      Author: marchi
 */

#ifndef ANALYZECHARGES_H_
#define ANALYZECHARGES_H_

#include <iostream>
#include "Metric.h"
#include "RhoAvg.h"
#include "CellSymmetry.h"
#include <cstdlib>
#include "time.h"


static Rho * ro=NULL;
static RhoDip * rodp0=NULL;
static int nframes=0;
static CellSymmetry * symmetry=NULL;
void compare(int nr,matrix co, matrix oc, rvec * x0, rvec * x1){
	rvec  x0a;
	rvec  x1a;
	rvec  xd;
	for(int i=0;i<nr;i++){
		x0a[XX]=oc[XX][XX]*x0[i][XX]+oc[XX][YY]*x0[i][YY]+oc[XX][ZZ]*x0[i][ZZ];
		x0a[YY]=oc[YY][XX]*x0[i][XX]+oc[YY][YY]*x0[i][YY]+oc[YY][ZZ]*x0[i][ZZ];
		x0a[ZZ]=oc[ZZ][XX]*x0[i][XX]+oc[ZZ][YY]*x0[i][YY]+oc[ZZ][ZZ]*x0[i][ZZ];

		x1a[XX]=oc[XX][XX]*x1[i][XX]+oc[XX][YY]*x1[i][YY]+oc[XX][ZZ]*x1[i][ZZ];
		x1a[YY]=oc[YY][XX]*x1[i][XX]+oc[YY][YY]*x1[i][YY]+oc[YY][ZZ]*x1[i][ZZ];
		x1a[ZZ]=oc[ZZ][XX]*x1[i][XX]+oc[ZZ][YY]*x1[i][YY]+oc[ZZ][ZZ]*x1[i][ZZ];
		xd[XX]=x0a[XX]-x1a[XX];
		xd[YY]=x0a[YY]-x1a[YY];
		xd[ZZ]=x0a[ZZ]-x1a[ZZ];
		xd[XX]=xd[XX]-rint(xd[XX]);
		xd[YY]=xd[YY]-rint(xd[YY]);
		xd[ZZ]=xd[ZZ]-rint(xd[ZZ]);
		cout << i << " " << xd[XX] << " " << xd[YY] << " " << xd[ZZ] << endl;
		if(i < 3){
			cout << i << " " << x0a[XX] << " " << x0a[YY] << " " << x0a[ZZ] << endl;
			cout << i << " " << x1a[XX] << " " << x1a[YY] << " " << x1a[ZZ] << endl;

		}
	}

}


static int AnalyzePolar(t_topology *top, t_trxframe *fr, t_pbc *pbc,
              int nr, gmx_ana_selection_t *sel[], void *data)
{
    /* Here, you can do whatever analysis your program requires for a frame. */
	Averages * Avgs=static_cast<Averages *> (data);

	RhoAvg<Rho> ** roavg0=Avgs->Ro;
    RhoAvg<Rho> * & roavg=*roavg0;
    RhoAvg<RhoDip> * rodipavg0=Avgs->RoDip;
    RhoAvg<RhoDip> & rodipavg=*rodipavg0;

    srand ( time(NULL) );

    Grid<1>::set(fr->box);
	Grid<DIM>::set(fr->box);

    if(!symmetry) symmetry=new CellSymmetry(top,fr,pbc);
    if(!ro) {
    	ro=new Rho [ngrps];
    }
    if(!rodp0) {
    	rodp0=new RhoDip;
    }
    RhoDip & rodp=*rodp0;
    Metric Met;
    Met(fr->box);
/*
    rvec * xx0=new rvec [natoms];
    rvec * xx1=new rvec [natoms];
*/

    if(bWater){
    	symmetry->Mypbc();

    	cidx[Tot].MolTot(MyMols);

/*
    	MyDip0.setCoord(Met,fr->x,cidx[Tot]);

    	for(int i=0;i<MyDip0.getNR();i++){
    		xx0[i][XX]=MyDip0[i][XX];
    		xx0[i][YY]=MyDip0[i][YY];
    		xx0[i][ZZ]=MyDip0[i][ZZ];
    	}
*/

/*
    	for(int o=0;o<200;o++){
    		int n= rand() % MyMols.nr;
*/
    	for(int n=0;n<200;n++){
    		cout << " Illa " << n << endl;
    		cidx[Slt].Mol(n,MyMols);
    		symmetry->Center1(cidx[Slt]);
    		cidx[Sol].MolCompl(n,MyMols);

    		chg[Slt].setCoord(Met,fr->x,cidx[Slt]);
    		chg[Sol].setCoord(Met,fr->x,cidx[Sol]);
    		chg[Tot].setCoord(Met,fr->x,cidx[Tot]);

    		MyDip.setCoord(Met,fr->x,cidx[Sol]);

/*
    		for(int i=0;i<MyDip.getNR();i++){
    			xx1[i][XX]=MyDip[i][XX];
    			xx1[i][YY]=MyDip[i][YY];
    			xx1[i][ZZ]=MyDip[i][ZZ];
    		}
    		compare(MyDip.getNR() ,Met.getCO(), Met.getOC(),xx0,xx1);
*/

/*
    		cout << " Ulla " << endl;
    		for(int i=0;i<MyDip.getNR();i++){
    			cout << i << " " << xx0[i][XX]-MyDip[i][XX] << " " << xx0[i][YY]-MyDip[i][YY] << " " <<  xx0[i][ZZ]-MyDip[i][ZZ] << endl;
    		}
    		cout << " Ulla " << endl;
*/

    		rodp.Density(MyDip);


    		ro[Slt].Density(chg[Slt]);
    		ro[Sol].Density(chg[Sol]);
    		ro[Tot]=ro[Slt]+ro[Sol];


    		roavg[Slt].Accumulate(ro[Slt]);
    		roavg[Sol].Accumulate(ro[Sol]);
    		roavg[Tot].Accumulate(ro[Tot]);
    		rodipavg.Accumulate(rodp);
    		nframes++;
    	}
    } else {
    	symmetry->Center();

//    symmetry->Fit();


    	for(int i=0;i<ngrps;i++)
    		chg[i].setCoord(Met,fr->x,cidx[i]);


    	MyDip.setCoord(Met,fr->x,cidx[Sol]);

    	rodp.Density(MyDip);

    	ro[Slt].Density(chg[Slt]);
    	ro[Sol].Density(chg[Sol]);
    	ro[Tot]=ro[Slt]+ro[Sol];

    	roavg[Slt].Accumulate(ro[Slt]);
    	roavg[Sol].Accumulate(ro[Sol]);
    	roavg[Tot].Accumulate(ro[Tot]);
    	rodipavg.Accumulate(rodp);
    	nframes++;
    }
    return 0;
}

#endif /* ANALYZECHARGES_H_ */
